MCMC convergence diagnostics
Adam Howes
Abstract
Background In the report naomi-simple_fit with parameter tmbstan = TRUE, we used the NUTS algorithm to perform MCMC inference for the simplified Naomi model.
Task Here we assess whether or not the results of the MCMC are suitable using a range of diagnostic tools.
We start by obtaining results from the latest version of naomi-simple_fit with tmbstan = TRUE.
out <- readRDS("depends/out.rds")
mcmc <- out$mcmc$stanfitThis MCMC took 3.28 days to run
cbpalette <- c("#56B4E9", "#009E73", "#E69F00", "#F0E442", "#0072B2", "#D55E00", "#CC79A7", "#999999")
bayesplot::color_scheme_set("viridisA")
ggplot2::theme_set(theme_minimal())1 \(\hat R\)
We are looking for values of \(\hat R\) less than 1.1 here.
rhats <- bayesplot::rhat(mcmc)
bayesplot::mcmc_rhat_data(rhats) %>%
ggplot(aes(x = value, y = parameter, col = description)) +
geom_segment(aes(yend = parameter, xend = ifelse(min(value) < 1, 1, -Inf)), na.rm = TRUE, alpha = 0.7) +
scale_color_manual(values = "#E69F00") +
geom_vline(xintercept = 1.05, linetype = "dashed", col = "grey40") +
labs(x = "Potential scale reduction factor", y = "NUTS parameter", col = "") +
theme_minimal() +
theme(
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
panel.grid.major = element_blank()
)
ggsave("rhat.png", h = 3, w = 6.25)
(big_rhats <- rhats[rhats > 1.1])## named numeric(0)length(big_rhats) / length(rhats)## [1] 02 ESS ratio
Reasonable to be worried about values less than 0.1 here.
ratios <- bayesplot::neff_ratio(mcmc)
breaks <- c(0, 0.1, 0.25, 0.5, 0.75, 1)
bayesplot::mcmc_neff_data(ratios) %>%
ggplot(mapping = aes(x = value, y = parameter, color = description)) +
geom_segment(aes(yend = parameter, xend = -Inf), na.rm = TRUE, alpha = 0.7) +
scale_color_manual(values = c("#56B4E9", "#009E73", "#E69F00")) +
geom_vline(xintercept = 0.1, linetype = "dashed", col = "grey40") +
geom_vline(xintercept = 0.5, linetype = "dashed", col = "grey40") +
geom_vline(xintercept = 1, linetype = "dashed", col = "grey40") +
labs(x = "ESS ratio", y = "NUTS parameter", col = "") +
theme_minimal() +
theme(
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
panel.grid.major = element_blank()
)
ggsave("ratio.png", h = 3, w = 6.25)3 ESS
What are the total effective sample sizes?
#' I think that this $summary should be all of the chains grouped together
mcmc_summary <- summary(mcmc)$summary
data.frame(mcmc_summary) %>%
tibble::rownames_to_column("param") %>%
ggplot(aes(x = n_eff)) +
geom_histogram(alpha = 0.8) +
labs(x = "ESS", y = "Count")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggsave("ess.png", h = 3, w = 6.25)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.4 Autocorrelation
How much autocorrelation is there in the chains?
bayesplot::mcmc_acf(mcmc, pars = vars(starts_with("beta")))
5 Univariate traceplots
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("beta")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("logit")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("log_sigma")))
5.1 Prevalence model
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_rho_x[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_rho_xs[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("us_rho_x[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("us_rho_xs[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_rho_a[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_rho_as[")))
5.2 ART model
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_alpha_x[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_alpha_xs[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("us_alpha_x[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("us_alpha_xs[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_alpha_a[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_alpha_as[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("u_alpha_xa[")))
5.3 Other
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("ui_lambda_x[")))
5.4 ANC model
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("ui_anc_rho_x[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("ui_anc_alpha_x[")))
bayesplot::mcmc_trace(mcmc, pars = vars(starts_with("log_or_gamma["))) #' N.B. these are from the ANC attendance model
6 Pairs plots
6.1 Paired PCA parameters
Variation between units can be explained by high correlation and high standard deviation, or low correlation and low standard deviation. Hence there is an unidentifiabiility here that leads to correlated posteriors.
bayesplot::mcmc_pairs(mcmc, pars = c("log_sigma_alpha_as", "logit_phi_alpha_as"), diag_fun = "hist", off_diag_fun = "hex")
bayesplot::mcmc_pairs(mcmc, pars = c("log_sigma_rho_as", "logit_phi_rho_as"), diag_fun = "hist", off_diag_fun = "hex")
bayesplot::mcmc_pairs(mcmc, pars = c("log_sigma_alpha_a", "logit_phi_alpha_a"), diag_fun = "hist", off_diag_fun = "hex")
(plot <- bayesplot::mcmc_pairs(mcmc, pars = c("log_sigma_rho_a", "logit_phi_rho_a"), diag_fun = "hist", off_diag_fun = "hex"))
ggsave("rho_a.png", plot, h = 4, w = 6.25)6.2 ART attendance
There is a prior suspicion (from Jeff, Tim, Rachel) that the ART attendance model is unidentifiable.
Let’s have a look at the pairs plot for neighbouring districts and the log_or_gamma parameter.
area_merged <- sf::read_sf(system.file("extdata/demo_areas.geojson", package = "naomi"))
nb <- area_merged %>%
filter(area_level == max(area_level)) %>%
bsae::sf_to_nb()
neighbours_pairs_plot <- function(par, i) {
neighbour_pars <- paste0(par, "[", c(i, nb[[i]]), "]")
bayesplot::mcmc_pairs(mcmc, pars = neighbour_pars, diag_fun = "hist", off_diag_fun = "hex")
}
# area_merged %>%
# filter(area_level == max(area_level)) %>%
# print(n = Inf)Here are Nkhata Bay and neighbours:
neighbours_pairs_plot("log_or_gamma", 5) 
And here are Blantyre and neighbours:
neighbours_pairs_plot("log_or_gamma", 26)
7 NUTS specific assessment
np <- bayesplot::nuts_params(mcmc)
saveRDS(np, "nuts-params.rds")Are there any divergent transitions?
np %>%
filter(Parameter == "divergent__") %>%
summarise(n_divergent = sum(Value))bayesplot::mcmc_nuts_divergence(np, bayesplot::log_posterior(mcmc))
We can also use energy plots (Betancourt 2017): ideally these two histograms would be the same When the histograms are quite different, it may suggest the chains are not fully exploring the tails of the target distribution.
bayesplot::mcmc_nuts_energy(np)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Original computing environment
sessionInfo()## R version 4.2.0 (2022-04-22)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS 13.3.1
##
## Matrix products: default
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] multi.utils_0.1.0 sf_1.0-9 bayesplot_1.9.0 patchwork_1.1.2 inf.utils_0.1.0 aghq_0.4.1 tmbstan_1.0.4 rstan_2.21.5 StanHeaders_2.21.0-7 TMB_1.9.2
## [11] Matrix_1.5-4.1 forcats_0.5.2 stringr_1.5.0 dplyr_1.0.10 purrr_1.0.1 readr_2.1.3 tidyr_1.3.0 tibble_3.2.1 ggplot2_3.4.0 tidyverse_1.3.1
##
## loaded via a namespace (and not attached):
## [1] colorspace_2.0-3 deldir_1.0-6 ellipsis_0.3.2 class_7.3-20 ggridges_0.5.3 rprojroot_2.0.3 fs_1.6.1 rstudioapi_0.14 proxy_0.4-27 hexbin_1.28.2
## [11] farver_2.1.1 fansi_1.0.4 mvtnorm_1.1-3 lubridate_1.8.0 xml2_1.3.3 codetools_0.2-18 splines_4.2.0 cachem_1.0.6 knitr_1.41 polynom_1.4-1
## [21] jsonlite_1.8.4 broom_0.8.0 naomi_2.8.5 dbplyr_2.1.1 compiler_4.2.0 httr_1.4.4 backports_1.4.1 assertthat_0.2.1 fastmap_1.1.0 cli_3.6.1
## [31] s2_1.1.1 htmltools_0.5.3 prettyunits_1.1.1 tools_4.2.0 gtable_0.3.1 glue_1.6.2 reshape2_1.4.4 posterior_1.2.2 wk_0.7.0 V8_4.2.2
## [41] Rcpp_1.0.10 jquerylib_0.1.4 cellranger_1.1.0 vctrs_0.6.1 spdep_1.2-7 mvQuad_1.0-6 tensorA_0.36.2 xfun_0.37 ps_1.7.3 rvest_1.0.2
## [51] bsae_0.2.7 lifecycle_1.0.3 statmod_1.4.36 orderly_1.4.3 ids_1.0.1 MASS_7.3-56 scales_1.2.1 ragg_1.2.2 hms_1.1.2 parallel_4.2.0
## [61] inline_0.3.19 yaml_2.3.7 curl_5.0.0 gridExtra_2.3 sass_0.4.4 loo_2.5.1 stringi_1.7.8 highr_0.9 checkmate_2.1.0 e1071_1.7-12
## [71] boot_1.3-28 pkgbuild_1.3.1 spData_2.2.1 rlang_1.1.0 pkgconfig_2.0.3 systemfonts_1.0.4 matrixStats_0.62.0 distributional_0.3.0 evaluate_0.20 lattice_0.20-45
## [81] labeling_0.4.2 cowplot_1.1.1 processx_3.8.0 tidyselect_1.2.0 traduire_0.0.6 bookdown_0.26 plyr_1.8.8 magrittr_2.0.3 R6_2.5.1 generics_0.1.3
## [91] DBI_1.1.3 pillar_1.9.0 haven_2.5.0 withr_2.5.0 units_0.8-0 sp_1.5-1 abind_1.4-5 modelr_0.1.8 crayon_1.5.2 uuid_1.1-0
## [101] KernSmooth_2.23-20 utf8_1.2.3 tzdb_0.3.0 rmarkdown_2.18 grid_4.2.0 readxl_1.4.0 isoband_0.2.6 data.table_1.14.6 callr_3.7.3 reprex_2.0.1
## [111] digest_0.6.31 classInt_0.4-8 numDeriv_2016.8-1.1 textshaping_0.3.6 openssl_2.0.5 RcppParallel_5.1.5 stats4_4.2.0 munsell_0.5.0 bslib_0.4.1 askpass_1.1Bibliography
Betancourt, Michael. 2017. “A Conceptual Introduction to Hamiltonian Monte Carlo.” arXiv Preprint arXiv:1701.02434.